feedback

Losing track of variable definitions

  • \(\alpha\) and \(\beta\) are confusing. Need to think about how to frame them better in the slides.
  • Options include:
    • footnotes/box in the corner with definitions for reminder
    • always writing them as the ratios
  • I wonder if it would be clearer to abandom \(S_1, S_0\) and instead spell out \(\Pr(\text{Symptomatic}\,|\,\text{untested})\), \(\Pr(\text{Asymptomatic}\,|\,\text{untested})\) (though clarify that’s not exactly what they mean)

Evan’s question on dependency of \(\alpha\) and \(\beta\)

  • \(\alpha = \dfrac{P(\text{test}_+| \text{untested}, S_1)}{P(\text{test}_+|\text{tested})}\)
  • \(\beta = \dfrac{P(\text{test}_+| \text{untested}, S_0)}{P(\text{test}_+|\text{tested})}\)
  • denominator is not treated as a random variable – this is the test positivity we observe

Notes on Independence (Pre-melding)

The number of events included here makes thinking about dependencies difficult, particularly since \(\alpha\) and \(\beta\) are related, but we sample them independently. It’s more clear to ignore the denominator since \(P(\text{test}_+|\text{tested})\) is just the observed test positivity and is a constant.

Observing \(P(\text{test}_+|\text{untested},S_0)\) doesn’t tell us anything directly about \(P(\text{test}_+|\text{untested},S_1)\). For example, if \(P(\text{test}_+|\text{untested},S_0)=1\), \(P(\text{test}_+|\text{untested},S_1)\) could be zero, or it could be 1, or anything in between.

If we had information about \(P(\text{test}_+|\text{untested}), P(S_1),\) and \(P(S_0)\), we could relate these quantities as

\[P(\text{test}_+|\text{untested}) = P(\text{test}_+, S_0|\text{untested}) P(S_0) + P(\text{test}_+, S_0|\text{untested})P(S_1).\]

We can also see their relationship with \(P(S_1|\text{untested}) = 1 - P(S_0|\text{untested})\) by considering

\[P(\text{test}_+|\text{untested}, S_0) = \dfrac{P(\text{test}_+, \text{untested}, S_0)}{P(\text{untested},S_0)} = \dfrac{P(\text{test}_+, S_0| \text{untested})}{P(S_0|\text{untested})},\] and similarly

\[P(\text{test}_+|\text{untested}, S_1) = \dfrac{P(\text{test}_+, S_1| \text{untested})}{P(S_1|\text{untested})} = \dfrac{P(\text{test}_+, S_1| \text{untested})}{1-P(S_0|\text{untested})},\]

but that said, with information on \(P(S_1|\text{untested})\) alone, we will would not be able to relate \(P(\text{test}_+|\text{untested}, S_0)\) and \(P(\text{test}_+|\text{untested}, S_1)\) directly because we wouldn’t know how the numerators \(P(\text{test}_+, S_0| \text{untested}), \;\; P(\text{test}_+, S_1| \text{untested})\) relate.

\[M(\theta) = \dfrac{\beta(1- P(S_1|\text{untested}))}{\beta(1-P(S_1|\text{untested})) + \alpha P(S_1|\text{untested})}\]

Summary: \(\alpha\) and \(\beta\) are dependent, but they are dependent on information we don’t have, so without this knowledge we sample them independently.

Notes on Independence (Post-melding)

During melding, we resample from \(\alpha, \beta, P(S_1|\text{untested})\) with the same weights, placing more weight on observations that yield observations with high density in the pooled distribution. If we compute \(M(\theta)\) with the constrained distributions of \(\theta\), the resulting distribution is the melded distribution.

When we sample from them from the melded priors for use in the probabilistic bias analysis, we account for this dependency by sampling them together.

source(here('analysis/base_functions/get_melded.R'))

source(here('analysis/base_functions/base_functions.R'))


prior_params <- list(
    alpha_mean = .95,
    alpha_sd = 0.08,
    alpha_bounds = NA,
   # alpha_bounds = c(.8,1),
    beta_mean = .15,
    beta_sd =.09,
    beta_bounds = NA,
  #  beta_bounds = c(0.002, 0.4),
    s_untested_mean = .03,
    s_untested_sd = .0225,
  #  s_untested_bounds = c(0.0018, Inf),
    s_untested_bounds = NA,
    p_s0_pos_mean = .4,
    # p_s0_pos_sd = .1225,
  p_s0_pos_sd = .1,
   p_s0_pos_bounds = NA,
  #  p_s0_pos_bounds = c(.25, .7),
    pre_nsamp = 1e6,
    post_nsamp = 1e5)

melded <- do.call(get_melded, prior_params)

melded$pre_melding %>%
  slice_sample(n=1e5) %>%
  mutate(source = "Before Melding") %>%
  bind_rows(melded$post_melding) %>%
  mutate(source = ifelse(is.na(source), "After Melding", source)) %>%
  ggplot(aes(x=alpha,y=beta,color = source)) +
  geom_point(alpha =.2,size=.5)+
    scale_color_manual(values = c(
      "Before Melding"= "#5670BF", 
      "Induced" = "#B28542",
      "After Melding" = "#418F6A")) +
  labs(color = "",
       x = TeX("$\\alpha$"),
       y = TeX("$\\beta$"),
       title = TeX("Relationship between $\\alpha$ and $\\beta$ Before and After Melding")) +
  theme_c(legend.position="top") +
  guides(color = guide_legend(override.aes = list(size = 5, alpha = 1)))

melded$post_melding %>%
  ggplot(aes(x=alpha,y=beta)) +
  geom_point(alpha =.2,size=.5)



melded <- do.call(get_melded, prior_params)


post <- melded$post_melding %>%
  select(P_A_testpos) %>%
  pivot_longer(everything())


melded$post_melding %>%
  mutate(induced_new = est_P_A_testpos(P_S_untested, alpha, beta)) %>%
  select(induced_new) %>%
  pivot_longer(everything()) %>%
  bind_rows(post) %>%
  ggplot(aes(x=value, fill =name)) +
  geom_density(alpha=.9)


melded$post_melding %>%
  mutate(induced_new = est_P_A_testpos(P_S_untested, alpha, beta)) %>%
  filter(induced_new != P_A_testpos)

melded_long <- reformat_melded(melded$post_melding,
                               melded$pre_melding,
                               pre_nsamp = prior_params$pre_nsamp,
                              p_s0_pos_mean = prior_params$p_s0_pos_mean,
                              p_s0_pos_sd = prior_params$p_s0_pos_sd,
                              p_s0_pos_bounds = NA)

################################
# FIRST WITHOUT MELDING 
################################
melded_long %>%
  filter(type == "Before Melding") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=FALSE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
  theme_c(legend.position = "top") +
  scale_fill_manual(values = list("Before Melding"= "#5670BF", "Induced" = "#B28542")) +
  labs(fill = "")

plt2 <-  melded_long  %>%
  filter(type!="After Melding" & name == "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=TRUE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
  theme_c(legend.position = "right",
          legend.text = element_text(size = 16)) +
  scale_fill_manual(values = list("Before Melding"= "#5670BF", "Induced" = "#B28542")) +
  labs(fill = "") 

cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))

ggsave(width=9, height = 8, "./presentation/figure/plot_no_melding.jpeg", dpi=300)


################################
# NOW WITH MELDING 
################################
plt1 <- melded_long %>%
  filter(name != "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=FALSE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
  theme_c(legend.position = "top") +
    scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  labs(fill = "")

plt2 <-  melded_long  %>%
  filter( name == "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=TRUE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
  theme_c(legend.position = "right",
          legend.text = element_text(size = 16)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  labs(fill = "") 

cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))

Evan’s question on why \(\alpha\) doesn’t change pre/post melding

Intuition: we form weights to sample from the log-pooled prior of \(\phi\), so since these weights aren’t heavily dependent on \(\alpha\), we are basically just randomly sampling from \(\alpha\).

prior_params <- list(
    alpha_mean = .95,
    alpha_sd = 0.08,
    alpha_bounds = NA,
   # alpha_bounds = c(.8,1),
    beta_mean = .15,
    beta_sd =.09,
    beta_bounds = NA,
  #  beta_bounds = c(0.002, 0.4),
    s_untested_mean = .03,
    s_untested_sd = .0225,
  #  s_untested_bounds = c(0.0018, Inf),
    s_untested_bounds = NA,
    p_s0_pos_mean = .4,
    # p_s0_pos_sd = .1225,
  p_s0_pos_sd = .1,
   p_s0_pos_bounds = NA,
  #  p_s0_pos_bounds = c(.25, .7),
    pre_nsamp = 1e6,
    post_nsamp = 1e5)


get_theta <- function(alpha_mean = 0.9,
                       alpha_sd = 0.04,
                       alpha_bounds = NA,
                       beta_mean = .15,
                       beta_sd =.09,
                       beta_bounds = NA,
                       s_untested_mean = .025,
                       s_untested_sd = .0225,
                       s_untested_bounds = NA,
                       p_s0_pos_mean = .4,
                       p_s0_pos_sd = .1225,
                       p_s0_pos_bounds = NA,
                       pre_nsamp = 1e4,
                      post_nsamp = 1e3,
                       include_corrected = TRUE,
                       bde = FALSE) {
  tibble(alpha = sample_gamma_density(pre_nsamp,
                                               mean = alpha_mean,
                                               sd = alpha_sd,
                                               bounds = alpha_bounds),
                  beta= sample_beta_density(pre_nsamp,
                                            mean = beta_mean,
                                            sd = beta_sd,
                                            bounds = beta_bounds),
                  P_S_untested = sample_beta_density(pre_nsamp,
                                                     mean = s_untested_mean,
                                                     sd = s_untested_sd,
                                                     bounds = s_untested_bounds)) %>%
    mutate(phi_induced = est_P_A_testpos(P_S_untested = P_S_untested,
                                         alpha = alpha,
                                         beta=beta))
  

}



tibble(alpha =seq(.7,1.3, length =100),
       beta = .15,
       P_S_untested = .03,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
  ggplot(aes(x=alpha, y = induced)) +
  geom_line() +
  theme_c() +
  labs(x= TeX("$\\alpha$"),
       y= TeX("$M(\\theta)$"),
       title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
       subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"))+
  ylim(.4, .95)



tibble(beta =seq(.05,.45, length =100),
       alpha = .95,
       P_S_untested = .03,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
  ggplot(aes(x=beta, y = induced)) +
  geom_line()+
  theme_c() +
  labs(x= TeX("$\\beta$"),
       y= TeX("$M(\\theta)$"),
        title = TeX("$M(\\theta)$ Across Values of $\\beta$"),
       subtitle = TeX("$P(S_1|untested),\\alpha$ are Fixed")) +
  ylim(.4, .95)



tibble(P_S_untested =seq(.05,.18, length =100),
       alpha = .95,
       beta=.15,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
  ggplot(aes(x=P_S_untested, y = induced)) +
  geom_line()+
  theme_c() +
  labs(x= TeX("$\\beta$"),
       y= TeX("$M(\\theta)$"),
        title = TeX("$M(\\theta)$ Across Values of $P(S_1|untested)$"),
       subtitle = TeX("$\\beta,\\alpha$ are Fixed")) +
  ylim(.4, .95)

tibble(alpha =seq(.7,1.3, length =100),
       beta = .15,
       P_S_untested = .03,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) %>%
  ggplot(aes(x=alpha, y = induced)) +
  geom_line() +
  theme_c() +
  labs(x= TeX("$\\alpha$"),
       y= TeX("$M(\\theta)$"),
       title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
       subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"))+
  ylim(.4, .95)

In Figure 1, we see increasing the mean of alpha corresponds to less density toward 1 and greater density at lower values of \(\phi\), but the effect is relatively small.

res <- map_df(seq(.8,1.3, by = .1), ~ {
  params <- prior_params
  params$alpha_mean <- .x
  theta <- do.call(get_theta, params)
  theta %>%
    mutate(alpha_mean = .x)
} )


res %>%
  mutate(alpha_mean = factor(alpha_mean)) %>%
  ggplot(aes(x=phi_induced, fill = alpha_mean)) +
  geom_density(alpha =.5) +
  viridis::scale_fill_viridis(option = "mako", discrete=TRUE, begin=.1, end = .9,direction=-1) +
  theme_c(legend.position="right",
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 14)) +
  labs(title = TeX("Induced Distribution on $\\phi$ by Mean of Alpha"),
       x = TeX("$\\phi^{induced}$"),
       fill = TeX("Mean of $\\alpha$")) +
  guides(fill = guide_legend(override.aes = list(alpha =1)))
\label{fig:alpha-mean}

Figure 1:

In particular, the effect of \(\alpha\) is limited by the small values of \(P(S_1|\text{untested})\). When \(P(S_1|untested)\) increases, \(\alpha\) makes a larger difference in the value of \(M(\theta)\) ( 2).

Red line corresponds to where \(P(S_1|\text{untested})\) is closest to its mean.

res <- map_df(seq(0.001, .2, length = 50), ~ {
  tibble(alpha =seq(.7,1.3, length =100),
       beta = .15,
       P_S_untested = .x,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) 
})

realistic <- res %>% 
  filter(abs(P_S_untested - 0.02536735) < 1e-4)


res %>%
  ggplot(aes(x=alpha, y = induced, color = P_S_untested, group =P_S_untested )) +
  geom_line() +
  geom_line(data = realistic,
            aes(x=alpha,
                y = induced,
                color = P_S_untested,
                group =P_S_untested ),
            color = 'red',
            linewidth =1.5) +
  labs(x= TeX("$\\alpha$"),
       y= TeX("$M(\\theta)$"),
       title = TeX("$M(\\theta)$ Across Values of $\\alpha$"),
       subtitle = TeX("$P(S_1|untested),\\beta$ are Fixed"),
       color = TeX("$P(S_1|untested)$"))+
  viridis::scale_color_viridis(option="magma", end=.9, direction = -1)+
  theme_c(legend.position="right",
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 14))



res <- map_df(seq(0.001, .2, length = 50), ~ {
  tibble(alpha =.95, 
       beta = seq(.03,.4, length =100),
       P_S_untested = .x,
       induced =  (beta * (1 - P_S_untested)) / (( beta * (1 - P_S_untested)) + (alpha * P_S_untested))) 
})


realistic <- res %>% 
  filter(abs(P_S_untested - 0.02536735) < 1e-4)


res %>%
  ggplot(aes(x=beta, y = induced, color = P_S_untested, group =P_S_untested )) +
  geom_line() +
  geom_line(data = realistic,
            aes(x=beta,
                y = induced,
                color = P_S_untested,
                group =P_S_untested ),
            color = 'red',
            linewidth =1.5) +
  labs(x= TeX("$\\beta$"),
       y= TeX("$M(\\theta)$"),
       title = TeX("$M(\\theta)$ Across Values of $\\beta$"),
       subtitle = TeX("$P(S_1|untested),\\alpha$ are Fixed"),
       color = TeX("$P(S_1|untested)$"))+
  viridis::scale_color_viridis(option="magma", end=.9, direction = -1)+
  theme_c(legend.position="right",
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 14))
\label{fig:effect}\label{fig:effect}

Figure 2:

Comparing Melded Distributions by Mean of \(\alpha\)

prior_params <- list(
    alpha_mean = .95,
    alpha_sd = 0.08,
    alpha_bounds = NA,
   # alpha_bounds = c(.8,1),
    beta_mean = .15,
    beta_sd =.09,
    beta_bounds = NA,
  #  beta_bounds = c(0.002, 0.4),
    s_untested_mean = .03,
    s_untested_sd = .0225,
  #  s_untested_bounds = c(0.0018, Inf),
    s_untested_bounds = NA,
    p_s0_pos_mean = .4,
    # p_s0_pos_sd = .1225,
  p_s0_pos_sd = .1,
   p_s0_pos_bounds = NA,
  #  p_s0_pos_bounds = c(.25, .7),
    pre_nsamp = 1e5,
    post_nsamp = 1e4)



res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
  
         melded <- do.call(get_melded, params)
       
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })


res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(\alpha\) SD

Increased SD

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$alpha_sd <- .3

         # params$s_untested_mean <- .2
         melded <- do.call(get_melded, params)
   
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })


res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Decreased SD

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$alpha_sd <- .04

         # params$s_untested_mean <- .2
         melded <- do.call(get_melded, params)
   
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })


res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(P(S_1|\text{untested})\) Mean

Upping \(P(S_1|\text{untested})\) to be .5

Unreasonable, but just to see

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$s_untested_mean <- .3

         # params$s_untested_mean <- .2
         melded <- do.call(get_melded, params)
   
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })


res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Upping \(P(S_1|\text{untested})\) to be .5

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$s_untested_mean <- .5

         # params$s_untested_mean <- .2
         melded <- do.call(get_melded, params)
   
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })


res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Comparing Melded Distributions by Mean of \(\alpha\) with Increased \(\alpha\) SD

Increased to .15

prior_params <- list(
    alpha_mean = .95,
    alpha_sd = 0.08,
  # alpha_sd = 0.15,
    alpha_bounds = NA,
   # alpha_bounds = c(.8,1),
    beta_mean = .15,
    beta_sd =.09,
    beta_bounds = NA,
  #  beta_bounds = c(0.002, 0.4),
    s_untested_mean = .03,
    s_untested_sd = .0225,
  #  s_untested_bounds = c(0.0018, Inf),
    s_untested_bounds = NA,
    p_s0_pos_mean = .4,
    # p_s0_pos_sd = .1225,
  p_s0_pos_sd = .1,
   p_s0_pos_bounds = NA,
  #  p_s0_pos_bounds = c(.25, .7),
    pre_nsamp = 1e5,
    post_nsamp = 1e4)




res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$alpha_sd <- .15

         params$s_untested_mean <- .03
         melded <- do.call(get_melded, params)
         # melded$post_melding %>%
         #   mutate(alpha_mean = .x) %>%
         #   cbind(alpha_before = melded$pre_melding$alpha)
         
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })



res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

# 
# res <- map_df(seq(.5, 2, length =5),
#       ~ {
#          params <- prior_params
#          params$alpha_mean <- .x
#          params$alpha_sd <- .15
# 
#          params$s_untested_mean <- .2
#          melded <- do.call(get_melded, params)
#          # melded$post_melding %>%
#          #   mutate(alpha_mean = .x) %>%
#          #   cbind(alpha_before = melded$pre_melding$alpha)
#          
#          reformat_melded(melded$post_melding,
#                          melded$pre_melding, 
#                          params$pre_nsamp,
#                           p_s0_pos_mean = params$p_s0_pos_mean,
#                             p_s0_pos_sd = params$p_s0_pos_sd,
#                             p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
#            mutate(alpha_mean = .x)
#          
#        })
# 
# 
# res %>%
#   mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
#  # filter(name == "$\\alpha$") %>%
#   ggplot(aes(x = value, fill = type)) +
#   geom_density(alpha=.8) +
#   facet_grid(name~alpha_mean,
#              scales = "free_y",
#              labeller=as_labeller(TeX,
#                                   default=label_parsed)) +
#   scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
#                                  "After Melding" = "#418F6A")) +
#   theme_c() +
#   labs(fill="")
#

Increased to .25

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$alpha_sd <- .25

         params$s_untested_mean <- .03
         melded <- do.call(get_melded, params)
         # melded$post_melding %>%
         #   mutate(alpha_mean = .x) %>%
         #   cbind(alpha_before = melded$pre_melding$alpha)
         
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })



res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Fixing \(\beta,P(S_1|\text{untested})\)

res <- map_df(seq(.5, 2, length =5),
      ~ {
         params <- prior_params
         params$alpha_mean <- .x
         params$beta_sd <- 0
         params$s_untested_sd <- 0
         melded <- do.call(get_melded, params)
         # melded$post_melding %>%
         #   mutate(alpha_mean = .x) %>%
         #   cbind(alpha_before = melded$pre_melding$alpha)
         
         reformat_melded(melded$post_melding,
                         melded$pre_melding, 
                         params$pre_nsamp,
                          p_s0_pos_mean = params$p_s0_pos_mean,
                            p_s0_pos_sd = params$p_s0_pos_sd,
                            p_s0_pos_bounds = params$p_s0_pos_bounds) %>%
           mutate(alpha_mean = .x)
         
       })



res %>%
  mutate(alpha_mean = factor(round(alpha_mean,1))) %>%
 # filter(name == "$\\alpha$") %>%
  ggplot(aes(x = value, fill = type)) +
  geom_density(alpha=.8) +
  facet_grid(name~alpha_mean,
             scales = "free_y",
             labeller=as_labeller(TeX,
                                  default=label_parsed)) +
  scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  theme_c() +
  labs(fill="")

Fixed versus Free Scales

Even allowing for different scales, \(\alpha\) does not change substantially.

Fixed Scale

melded <- do.call(get_melded, prior_params)


post <- melded$post_melding %>%
  select(P_A_testpos) %>%
  pivot_longer(everything())


melded_long <- reformat_melded(melded$post_melding,
                               melded$pre_melding,
                               pre_nsamp = prior_params$pre_nsamp,
                              p_s0_pos_mean = prior_params$p_s0_pos_mean,
                              p_s0_pos_sd = prior_params$p_s0_pos_sd,
                              p_s0_pos_bounds = NA)


################################
# NOW WITH MELDING 
################################
plt1 <- melded_long %>%
  filter(name != "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=FALSE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=3) +
  theme_c(legend.position = "top") +
    scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  labs(fill = "")

plt2 <-  melded_long  %>%
  filter( name == "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=TRUE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4) +
  theme_c(legend.position = "right",
          legend.text = element_text(size = 16)) +
scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  labs(fill = "") 

cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))

Free Scale

################################
# NOW WITH MELDING 
################################
plt1 <- melded_long %>%
  filter(name != "$P(S_0|test+,untested)$") %>%
  ggplot(aes(x=value, fill = type)) +
  geom_density( alpha = .8, show.legend=FALSE) +
  facet_wrap(~name, labeller=as_labeller(TeX, default = label_parsed), ncol=4, scales="free_y") +
  theme_c(legend.position = "top") +
    scale_fill_manual(values = c("Before Melding"= "#5670BF", "Induced" = "#B28542",
                                 "After Melding" = "#418F6A")) +
  labs(fill = "")


cowplot::plot_grid(plt1,plt2, ncol=1, rel_widths = c(1,.9))

Conclusions

The lack of an effect of \(\alpha\) is primarily a result of the functional form of \(M\); because \(\beta\) has a larger effect

M <- function(input1, input2) {
 # x <- 0.5
  x <- 0.03
  input1*(1-x) / (input1*(1-x) + input2*x)
}

M  <- Vectorize(M )

inp1 <- seq(0.01,.4, length = 50) 
names(inp1) <- inp1
inp2 <- seq(.7,1.3, length = 50)  
names(inp2) <- inp2
z <- outer(inp1,inp2, FUN =M)


# 3D plot
par(mar = c(0, 0, 0, 0)) 
persp(inp1, inp2, z,  xlab = "beta", ylab= "alpha", theta = 60, phi = 30)

library(plotly)



ax <- list(
  title = "alpha"
)


ay <- list(
  title = "beta"
)




fig <- plot_ly(z=~z)

fig <- plot_ly( x= ~colnames(z), 
                y =~ rownames(z),  
                z=~z)

fig %>% add_surface() %>% layout(scene = list(xaxis = ax,
                                              yaxis = ay))

Aaron’s question on notation

Evan’s question on \(P(S_0|\text{test}_+, \text{untested}) \neq P(S_0|\text{test}_+, \text{tested})\)

The fact that \(P(S_0|\text{test}_+, \text{untested}) > P(S_0|\text{test}_+, \text{tested})\) is something important to think about. It has been a while since my initial search for meta-analyses, so another look would be good to get a better sense for the estimates from screening studies and the extent to which they differ.